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ABSTRACT 



Planet formation occurs within the gas- and dust-rich environments of protoplane- 
tary disks. Observations of these objects show that the growth of primordial submicron- 
sized particles into larger aggregates occurs at the earliest evolutionary stages of the 
disks. However, theoretical models of particle growth that use the Smoluchowski equa- 
tion to describe collisional coagulation and fragmentation have so far failed to produce 
large particles while maintaining a significant population of small grains. This has been 
generally attributed to the existence of two barriers impeding growth due to bouncing 
and fragmentation of colliding particles. In this paper, we demonstrate that the impor- 
tance of these barriers has been artificially inflated through the use of simplified models 
that do not take into account the stochastic nature of the particle motions within the 
gas disk. We present a new approach in which the relative velocities between two par- 
ticles is described by a probability distribution function that models both deterministic 
motion (from the vertical settling, radial drift and azimuthal drift) and stochastic mo- 
tion (from Brownian motion and turbulence). Taking both into account can give quite 
different results to what has been considered recently in other studies. We demonstrate 
the vital effect of two "ingredients" for particle growth: the proper implementation of a 
velocity distribution function that overcomes the bouncing barrier and, in combination 
with mass transfer in high-mass-ratio collisions, boosts the growth of larger particles 
beyond the fragmentation barrier. A robust result of our simulations is the emergence 
of two particle populations (small and large), potentially explaining simultaneously a 
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number of long-standing problems in protoplanetary disks, including planetesimal for- 
mation close to the central star, the presence of mm- to cm-size particles far out in the 
disk, and the persistence of micron-size grains for millions of years. 

Subject headings: accretion, accretion disks - solar system: formation 



1. Introduction 



1.1. Observations confront theory 

Protoplanetary disks are universally thought to be the birthplace of planets. Furthermore, the 



ubiquity of terrestrial to gas-giant exoplanets discovered around other stars (Cassan et al. (2012); 



Howard et al. (2012), see also the compilation at http://www.circumstellardisks.org) strongly sug- 



gests that the planet formation process must be robust and efficient. This simple fact, however, 
has been until now rather difficult to reconcile with a number of observational properties of proto- 
planetary disks. 

Consider for instance a scenario in which planets are created by the gradual coagulation of small 
dust grains into progressively larger particles, eventually leading to the formation of planetesimals 
and then planets. This idea turns out to be surprisingly hard to model within the constraints 
set by disk observations. Indeed, signatures of ongoing gas accretion onto the host star typically 



disappear after lOMyr (Calvet et al. 2000). In the core-accretion scenario (Pollack et al. 1996), 



coagulation must be fast enough to drive growth from submicron-sized dust grains to earth-sized 
rocky planets while gas is still present in sufficient quantities to form gas giants. 

More crucially, particles as they grow must necessarily pass through a critical phase where their 
radial velocity (induced by gas drag) approaches a small but significant fraction of their orbital 
speed ( Weidenschilling| 1977 ) . They are prone to drift into the central star within 10 to a 100 orbits 
unless growth through that phase is rapid enough. In short, coagulation must be very efficient for 
planets to exist. 

However, protoplanetary disks are usually detected via a characteristic near- and mid-infrared 



emission that is in excess of the expected pure photospheric flux of their host stars ( Cohen & Kuhi 



1979). Objects showing such a feature in their spectral energy distribution are known as low- mass 



T-Tauri or intermediate-mass Herbig AeBe stars (Gillett & Stein 1971; Strom et al. 1972). The 



main contribution to the excess emission comes from hot micron-sized dust grains located within a 
few AU of the star. The coagulation efficiency of these tiny dust grains, upon collision, is usually 



thought to be much higher than that of larger bodies ( Dullemond Sz Dominik 2005 ) . This implies 



that the growth process would quickly render them invisible at these wavelengths, shifting the 



excess emission to longer wavelengths, in contrast to what is observed (Haisch et al. 2001 Cieza 
eTaI][2007] |Su et al.|[2009| ). 
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A very similar problem was raised with the observational discovery of "large" grains at large 
orbital radii away from the central star - more precisely, mm-to-cm-size grains, at radii at least as 
far out as a few tens of AU ( Holland et al.||1998 Wilner et al.||2005 Rodmann et al.||2006 Lommen 
et al,|[2009 ). Forming such grains at these orbital distances, where the dynamical timescale is more 
than a hundred times longer than around 1AU, is just as difficult as forming meter-sized objects in 
the inner stellar system. Furthermore, mm- and cm-size grains at 30AU are also subject to rapid 
inward drift. Again, coagulation has to be very efficient for growth to cm-size to happen despite 
these two problems, but cannot a priori be too efficient otherwise the outer (and inner) disk would 



disappear entirely on a much shorter timescale than the one observed (Garaud 2007; Hernandez 



et al. 2008). 



A possible way out of this impasse was proposed by Dullemond & Dominik (2005), who in- 



troduced the possibility of fragmentation in addition to coagulation. Fragmentation is a common 
outcome of high- velocity collisions (see review by Blum Wurm|[2008 ) . As long as the fragmenta- 
tion rate is high enough, and the size distribution of the fragments is sufficiently steep to favor the 
creation of small grains upon impact, then it is possible to maintain a healthy micron-size grain 
population at all times in the disk ( Dullemond Dominik|2005 Brauer et al.|2008 ). Unfortunately, 
while fragmentation solves one problem it re-introduces another - particle growth beyond a new 
"fragmentation barrier", corresponding to cm-size at a few AU and down to less than mm-size at 
30 AU, seems to be impossible unless additional physics are taken into account ( Brauer et al.j [2008 



Johansen et al. 2008). 



1.2. Modeling particle growth using ensemble models 

Upon realization of the true complexity of this problem, much effort was dedicated to modeling 
the collisional dynamics of particles in protoplanetary disks in more detail. The ensemble approach 
traditionally used to study aerosol growth in planetary atmospheres ( Kovetz &z Qlund|1969 Podolak 



& Podolak 1980) and grain growth in molecular clouds (Rossi et al. 1991 Ossenkopf 1993), in 



which one models the evolution of the particle size distribution function with the Smoluchowski 



coagulation/fragmentation equations (Smoluchowski 1916 Melzak 1957), has become one of the 
most commonly used techniques applied to the problem. We shall adopt it here too. 

The Smoluchowski equations are integro-differential equations. Information on the collisional 
dynamics of particles must be provided in the form of "collisional kernels" , which summarize all 
aspects of the collisions, including the relative velocities of the colliding particles, their cross sec- 
tions, their sticking and/or fragmentation probabilities, the fragment distribution, etc... Formally 
speaking, the coagulation kernel of two particles and "j" (the latter being two indices that 
reference the particle masses mi and rrij for instance) should be written as: 



K, 



dSi 



dSj I dl 



dAi 



■jHj] 
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p e (ef 3 -|Si,Sj, Aij,<5y^D)pA(Ai3|Sj,SjsD)p/(I|Si,Sj;D)ps.(Si;D)p Sj .(Sj;D)} (1) 
or, in other words, a high- dimension multivariate integral over 

• all possible variations in each particle's internal properties that result in the same parti- 
cle mass, loosely summarized here as the vector S, which may include shape, composition, 
porosity, or other, and is drawn from a multivariate probability distribution function (p.d.f. 
hereafter) pg. Note that if the masses of the two particles are different then the p.d.f.s for 
each of the particles may also be different. Since the grain composition could depend on 
the local temperature, the p.d.f.s should also depend on the local disk properties, loosely 
described here as the vector D. Finally, pg may also be a function of time, if the repeated 
collisions lead to internal changes in the particles - this has not be explicitly written here for 
simplicity. 

• all possible impact configurations, loosely summarized here as the vector I, which includes 
impact parameter and solid angle, and drawn from a p.d.f. pj. If collisional outcomes are 
independent of I, then the integral separates and yields the mean area cross section of the 
collision ay. 

• all possible relative velocities Ay (drawn from a p.d.f. pa) for the collisions of the two 
particles. Note that since particle motion is induced by collisional/frictional interaction with 
the gas, this p.d.f. depend on the respective structure of both particles, and on the disk model 
considered. 

• all possible sticking efficiencies efj for the collisions of the two particles, drawn from a p.d.f. 
p e . This p.d.f depends on everything else, including the local disk properties. 

The integrand in the square brackets is the product of the collisional relative velocity Ay , and the 
sticking efficiency efj, discussed above. 

Clearly, equation ([TJ is far too complicated to be of practical use. Tanaka et al. (2005) and 



Dullemond & Dominik (2005) dramatically simplified the kernel expression to a product of three 



means, which may themselves only depend on mean quantities (and on the disk model): 

Ktj = Sy(Si, Sj-; D)Ay(S;, S i; D)e*-(Si, S i; Ay; D) , (2) 

where ay is the mean collisional cross section of a pair of particles of mass rm and rrij, assuming 
mean internal properties and S^; Ay is a mean relative velocity and ey is the mean sticking 
efficiency of the collisions. A similar equation is used to approximate the fragmentation kernel. 



1.3. Improving collisional dynamics models 



Having generally adopted this much simpler formulation, studies of the evolution of the particle 
size distribution function then focused on refining parametric prescriptions for the mean quantities 
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involved. Progress was made on all three fronts - in characterizing the collisional velocities, the 
particle structure (which determines the collisional cross-section) and collisional outcomes. 

Thanks to the vast advances in supercomputing, it is now possible to study the dynamics of 
a vast number of inertial particles interacting with the disk turbulence. Using this information, 



one can then infer statistical properties of individual particle velocities (Nelson & Gressel 2010 



Carballido et al.||2011[ ), and of pairwise relative velocities ( |Carballido et al.||2008[ |2010[ ). The latter 
have helped begin to validate a theoretical prescription for the rms relative particle velocity of 



colliding particles proposed by Ormel Sz Cuzzi (2007), which is commonly used in dust coagulation 
models to construct the quantity Aj,- used in equation pj). However, the study of the p.d.f.s of 



collisional velocities is still very much in its infancy (see Carballido et al. (2010) for a first attempt 
at the problem). 

In parallel, much work has also been done to improve our understanding of the outcome of 
particle collisions, and typical particle structures (the latter often depending on the collisional 
history) using both experimental and numerical methods. It is generally accepted that collisions 
involving only micron-size particles have very high sticking efficiency and that successive collisions 



result in the formation of fractal aggregates (see Ossenkopf 1993 also see the review by Blum & 



Wurm 2008). The collisions of these aggregates with one another (or with monomers) then leads 



to further sticking (with very high probability), but also to compaction (Dominik & Tielens 1997 



Blum & Wurm [2000 ) . The resulting objects are roughly mm- or sub- mm size grains that have 



varying degrees of porosity (Blum et al. 2006). In short, for any collision involving only sub-mm 
particles, one can usually adopt a mean sticking efficiency efj ~ 1; the calculation of the area 
cross-section as a function of the particle mass, by contrast, is more complicated if one wishes to 



take into account all aforementioned effects (Ossenkopf 1993). 



Collisions involving larger particles have much more varied outcomes, depending on their in- 
ternal properties, relative velocities, and mass ratio. Given the sheer dimensionality of parameter 
space, and the 35 orders of magnitude mass-range between mm-size grains and Earth-size objects, 
a comprehensive study of the problem is impossible. It is generally agreed that sticking continues 
to occur for low-velocity (alternatively low-energy) collisions, that intermediate velocity (energy) 
collisions may lead to bouncing with compaction ( |Blum Miinch|l993 Weidling et al.|[200 9), and 
that the outcome of high- velocity collisions very much depends on the mass ratio of the two par- 
ticles - with equal mass collisions being much more likely to lead to complete fragmentation than 
unequal mass collisions (see the review by Blum Wurm|2008 ). In fact, high- velocity collisions in 



high-mass ratio events can also lead to partial or complete sticking (Teiser & Wurm 2009 Kothe 



et al. 2010). If the fraction of the projectile mass that sticks to the target is larger than the total 
amount of mass ejected away on impact - a phenomenon called "mass transfer" hereafter 



then 



high-mass-ratio collisions can lead to the growth of the target. 



However, the details of the transitions between each of these regimes is highly dependent on 
the particle masses and internal structure, on the impact parameter and on the collisional energy. 



- 6 - 



Guttler et al. (2010) recently proposed a complex model for sticking, bouncing and fragmentation 
outcomes that summarizes a suite of laboratory experiments. Numerical simulations of larger, 
centimeter-sized, dust aggregates show that the fragmentation velocities may be larger than the 
previously assumed lOOcm/s ( Geretshauser et aL||2011a ) and that the outcome may be dependent 
on the inhomogeneity of the aggregate (Geretshauser et al. 2011b), the target and projectile sizes 
as well as the porosity of the aggregates (Meru et al, in prep). 



1.4. Proposed work 

A number of papers have been published since 2005 which gradually take into account more 



and more of the aforementioned effects, as well as additional disk physics (Birnstiel et al. 2010), 
whilst keeping the kernel formulation given in equation 



Without giving an exhaustive list, we note in particular the work of Brauer et al. (2008) 



and Windmark et al. (2012a). Brauer et al. (2008) incorporate the collisional dust dynamics 



with a whole-disk model to study simultaneously the evolution of the particle size distribution 
function through coagulation and fragmentation as well as mean motions. They model the mean 
sticking and fragmentation probabilities using piecewise linear or piecewise parabolic functions, and 
consider sticking and fragmentation, but no bouncing. They also include effects such as "cratering" 
(or equivalently, erosioij^]), in which high-mass-ratio collisions only lead to the partial rather than 



the complete fragmentation of the larger body. Later on, Windmark et al. (2012a) discussed even 



more complex models, which include sticking, fragmentation, cratering and bouncing, taking into 



account the mass dependence of the regime boundaries (Weidling et al. 2012). Furthermore, they 



also considered the possibility of "mass transfer", where high-mass-ratio collisions lead to the net 



growth of the target. Windmark et al. (2012a) showed that the latter could lead to the growth 
of particles well-beyond the fragmentation barrier, but only if they are introduced - somewhat 
artificially - beyond it to start with. 

Nevertheless, despite all these added physics, most models fail to robustly produce planetesi- 
mals at 1AU and cm-size particles far out in the disk, for reasonable disk assumptions, whilst keeping 
a micron-size dust grain population consistent with observations. By contrast, the manner in which 
the coagulation and fragmentation kernels are constructed has not really been revisited since 2005. 
However, the statement that they can be approximated using products of three separately-taken 
averages is only technically correct if the quantities considered are mutually independent - which is 
certainly not the case. We demonstrate in this paper that the failure of previous models to reconcile 
disk observations with exoplanetary observations simply stems from the use of equation ^ , which 
vastly oversimplifies the problem and results in misleading conclusions concerning particle growth 



1 Note that 



Guttler et al. 



( 2010 1 refer to this mechanism as "erosion" while 



Windmark et al. 



1 2012a) use these 



terms interchangeably. 
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in protoplanetary disks. 



Our goal here is not to propose a complete new model for the construction of collisional kernels - 
certainly, none as complicated as the one suggested in equation ([I]). Instead, we aim to demonstrate 
that adding even just one element of stochasticity to the model, taking into account the full p.d.f. 
of relative velocities instead of their mean only, resolves a number of the aforementioned problems, 
at least qualitatively speaking. 



To our knowledge, Okuzumi et al. (2011) were the first to take into account velocity p.d.f.s 



explicitly in collisional dust growth models with application to protoplanetary disks. They focused 
on studying very small particles only, which interact via Brownian motion and mean drift. Here, 
we also consider larger particles, and include the effect of turbulence in driving stochastic motions. 
We loosely follow the same general methodology introduced in the ISIMA (International Summer 



Institute for Modeling in Astrophysics) report by Galvagni et al. (2011) and recently followed by 



Windmark et al. (2012b). We construct a plausible p.d.f. for the relative velocities of two particles 



interacting with the gaseous disk. We then use the latter to construct new collisional kernels, and 
study the resulting collisional particle growth and fragmentation. Our work nevertheless differs 
notably from previously published ideas, in several ways. First, we explicitly take into account 



the effects of turbulence in constructing the velocity p.d.f, by contrast with Okuzumi et al. (2011 ) 



Secondly, while Galvagni et al. (2011) and Windmark et al. (2012b) did take turbulence into 



account, they modeled the relative velocities p.d.f. as a Maxwellian, which neglects the difference 
between chaotically-driven particle motion (e.g. Brownian motion and turbulence) and regular 
particle motion (radial drift, azimuthal drift and vertical settling induced by gas drag). Moreover, 
the model proposed here does take drift into account more carefully, and is much more generally 
applicable. Finally, we point out and correct a mathematical error in the work of | Windmark et al.| 
( |2012b ) (which affects the construction of the mean sticking and fragmentation probabilities in the 
kernels). We provide here the correct derivation (within the assumptions made). 

In what follows, we present in Section [2] the general framework associated with the use and 
solution of the Smoluchowski equations. The construction of the kernels is discussed in Section |3j 
and contrasts the traditional approach with our new proposed one. Section [4] presents a sample disk 
model, and Section [5] shows how the evolution of the particle size distribution function in this disk 
dramatically changes between the two approaches. The principal finding is that particle growth 
beyond the traditional bouncing and fragmentation barriers is possible, and is a robust feature of all 
models that include relative velocity p.d.f.s in the kernel construction and consider the possibility 
of mass transfer in high-mass-ratio collisions. The particle size distribution function rapidly evolves 
into two populations, one composed of small dust grains and one composed of much larger particles. 
The quantitative details of these two populations, however, are strongly model-dependent. This is 
discussed in Section [6| together with recommendations for future work directions. 
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2. General framework 

We consider a local region of a protoplanetary disk, centered around the mid-plane at orbital 
radius r. The properties of the gas disk are assumed to be known, steady, and independent of the 
dust properties. This is a good approximation as long as the dust-to-gas mass ratio is not too large 
(which we will assume here), and as long as the evolution timescale of the dust size distribution 
function is short compared with the disk evolution timescale (which needs to be verified a posteriori 
for the disk model selected). 

The particle size (or mass) distribution varies with time, in a manner that is controlled both by 
the flux of particles in and out of the region, and through collisional coagulation or fragmentation. 



In this paper we neglect the former, for simplicity (the same approximation was used by Windmark 



et al. 2012b). We discuss the implications of this assumption in Section |6| In this local model, 
the evolution of the particle size distribution function is therefore simply governed by the Smolu- 
chowski coagulation-fragmentation equations. We now present these equations for completeness, 



but essentially follow the work of Brauer et al. (2008) 



2.1. Coagulation- Fragmentation equations 

Assuming that the masses of the dust particles can take any values among a discrete "mass- 
mesh" {rai,m2,...ra/}, we define iVj as the total number of particles of mass rrii. The coupled 
nonlinear evolution equations for the functions Ni(t) are the discrete form of the Smoluchowski 
coagulation- fragmentation equations ( Smoluchowski|[l916 Melzak 1957), in which, for i = 1..J, 



= \ E C^K^Nk - ^2 K hi N i N i 

j,k=l j=l 



1 

+ 2 E F jkNjN k N{ jk - £ FyNiNj . (3) 
j,k=l j=l 

In each of these equations, the first two terms on the right-hand-side model coagulation, while 
the last two terms model fragmentation. The quantities K^j and Fij are the coagulation and 
fragmentation kernels respectively, which contain crucial information on the relative velocities of 
particles of mass rrii and rrij, their collisional cross section, and their probability of sticking and 
fragmenting as discussed in Section [T] The construction of these kernels is discussed in Sections 2.2 
andO 

The first term on the right-hand-side of equation ^ models the increase in the number of 
particles of mass rrii resulting from the coagulation of two particles of mass rrij and mt such that 
rrij + rrik — mi. If the mass mesh is linearly spaced, then there exists a mass-point rrii which is 
exactly equal to rrij + m^, namely the point for which i = j + k. In that case, this coagulation term 
reduces to the more commonly-used g Yl % j=i Kj,i-jNjNi-j. 
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In order to follow the growth of particles across many orders of magnitude in mass, how- 
ever, one cannot realistically use a linearly-spaced mass-mesh - it is much more common to use a 
logarithmically-spaced one instead. With a non-uniform mesh, however, rrij + rrik usually falls in 
between two existing mass-points, with rrij_ < rrij + rrik < mj+. This leads to the introduction of 
the third-rank tensor Cjki, which assigns different probabilities for the coagulation of rrij and rrik 
into either mj_ or mj+. The details of how to construct Cjki in such a way as to conserve mass and 



coagulation rates were first discussed by |Kovetz & 01und| (|1969|) , and are very nicely summarized 
in 



Brauer et al. (2008). Following their work, we take 



Cjki = Pjk if i = i - 
Cjki = (1 - Pjk) if i = i + 
Cjki = otherwise 



where 



Pjk 



m. 



(rrij + m k ) 



rrii 



(4) 



The second and fourth term on the right-hand-side of equation ^ describe how particles of 
mass rrii disappear when they collide with a particle of mass rrij , and either coagulate (second term) 
or fragment (fourth term). 

Finally, the third term on the right-hand-side of equation ( 3 ) describes how N^ k particles of 
mass rm are created when two particles rrij and rrik (with rrij + rrik > rrii) collide and fragment. 
In what follows, we assume that in any fragmentation event the number of fragments generated is 
a power law with fixed index — £, where £ = 11/6 = 1.83 as in Brauer et al. (2008) for instance. 
Since the mass bins are logarithmically distributed, we have 



*&k = A jk rn]-tH(L 



(5) 



where H(L — i) is a discrete Heaviside function (defined so that H(L — i) = if i > L), L is the index 
of the mass-point mi, which is the mass of the largest fragment generated, and the normalization 
Ajk is uniquely defined by conservation of mass during the collision: 



rrij +m k = A jk ^ ■ 



2-f 



(6) 



We assume here for simplicity that mi is either equal to the largest mass-point still satisfying 
rriL < rrij + rrik, OT , i n the occasional very high- mass collisions events with rrij + rrik ^ m i (where 
mi is the largest mass-point in our mesh), mi is chosen to be equal to mi. 

It can be shown that this set of equations conserves the total mass of particles in the system 
exactly, with 

d 



i 

di^2 

i=l 



rrnNi = 



(7) 
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The coagulation-fragmentation equations are intrinsically nonlinear, and - aside from very 
specific kernels and initial conditions - must be evolved forward in time numerically. This system 
of equations is also inherently very stiff, so it is crucial to use an implicit time-stepping method. 
Here, we have chosen to use a simple Euler first-order time-centered scheme, and calculate the 
Jacobian of the right-hand-side explicitly. We use an automatic adaptive time-stepping scheme to 
ensure the required level of precision. The algorithm has been tested against well-known analytical 



solutions of these equations, reviewed by Aldous (1999). 



We continuously check for conservation of the total mass. As in Brauer et al. (2008), when the 
ratio of the smallest to the largest particle masses considered falls below the numerical precision 
of the system (which is the case in some of the simulations presented below), mass conservation 
becomes difficult to enforce. When the total mass drifts by more than one part in a million away 
from the initial mass, the difference is added back to the mass distribution function in the form 
of monomers at the smallest size considered. Comparison of the results using this simple fix, with 
those obtained using a quad-precision real arithmetic (which yields reliable results up to much 
larger particle sizes) shows that the difference is minimal. Since quad-precision real arithmetic 
requires about 10 times longer integration times, we use the fix in all simulations shown below. 

Finally, let us first discuss a subtle but rather important point concerning the discrete version 
of the Smoluchowski equations. In the derivations presented so far, we have defined Ni to be 
the number of particles of mass mj. This definition is useful for numerical purposes, but can be 
troublesome when interpreting the results. Indeed, consider for instance a uniformly distributed 
mass- mesh with a total of 100 mass-points, and a total of iVo particles. If the particle mass 
distribution function is itself uniform, then Ni = iVo/100 for all i. By contrast, had we selected 
200 mass-points instead, i.e. a finer mass-mesh, then Ni = Nq/200 for all i. This simple example 
illustrates that Ni is a resolution-dependent quantity. 

A much better approach is to consider the continuous particle mass distribution function 
dn/dm, defined so that, for any two masses m\ and m-2, f™ 2 (dn/dm)dm is the number of particles 
of mass between mi and rri2- This definition is now entirely independent of the mass-discretization 
used. In order to relate dn/dm to the numerically determined A^'s, simply note that, as long as 
there are many mass-points between mi and m,2, f™^(dn/dm)dm ~ X^i=»i ^» wnere i\ is the index 
of the closest mass-point to mi, and similarly for ii- Looking at the integral as a Riemann sum, 
one can then straightforwardly identify 



N — ^ n 
1 dm 



dn 

ami = ~j— 
dm 



(mi - mi_i) . (8) 

m=m,i 



In Section [5j we present all our results in terms of dn/dm; the latter are calculated from the 
discretized N using equation Q. 
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2.2. Ingredients for coagulation and fragmentation 

While the solution of the Schmoluchowski coagulation- fragmentation equation poses interesting 
mathematical and computational problems, most of the physical complexity of the system studied 
manifests itself in the construction of the coagulation and fragmentation kernels. 

The coagulation and fragmentation kernels Kij and Fij must be constructed for each (i, j) pair 
(i.e. for all collisions involving one particle of mass rrii and one particle of mass rrij). As discussed 
in Section [TJ they have traditionally been simplified into 

— — £ 

Kij = aijAijefj and Fij = Cy-Ay-e^ , (9) 

where ciij is the mean collisional cross-section of particle pairs (averaged over all possible 
individual internal structures resulting in the same mass) , A ^ is an estimate of the mean relative 
velocity of all pairs (averaged over all possible realizations), and ef^ is the mean sticking 

probability of a collisional event, given this mean structure and at that mean relative velocity (and 
similarly for e{j). This formulation is equivalent to assuming that all pairs have the same 
collisional cross-section (neglecting the possibility of variations in density through compaction, and 
variations in shape), collide with the same relative velocity A™ (ignoring the chaotic nature of 
particle motions), and that all these collisions have the same sticking and fragmentation probabil- 
ities efj and e\- (neglecting variations with impact parameter and inherent variability in collisional 
outcomes). 

Often used with specific oversimplified prescriptions for the fragmentation and sticking prob- 
abilities ef-(Ajj) and e^-(Ay), this formulation has, over the years, led to the introduction and 
discussion of arguably artificial "bouncing" and "fragmentation" barriers. As we shall demonstrate, 
however, taking into account even a single stochastic component - namely the chaotic nature of the 
particle motion - in the description of the collisions overcomes these barriers, and can explain the 
growth of large particles whilst retaining a population of small particles quite naturally. In what 
follows, we therefore still adopt a simple collisional cross-section model and a simple model for the 
fragmentation and sticking probabilities in individual collisions, but look at the effects of modeling 
the stochastic nature of the particle's relative velocities. We now describe these ingredients in more 
detail. 



Collisional cross-section. The collisional cross-section of two particles % and j is typically taken 
to be their area cross-section, although electrostatic or gravitational effects could be important for 
tiny dust particles and planetesimals respectively. The calculation of an area cross-section for two 
given particles requires knowledge of their shape and internal density, which both depend on their 
collisional history. Taking these effects into account properly can only be done via Monte Carlo 
simulations, as in Zsom & Dullemond (2008) for instance. In ensemble evolution models (such as 
the one presented here) one is forced instead to make simplifying assumptions. Following most 
previous work on the subject (Tanaka et al. 2005 Dullemond & Dominik 2005) we will consider 
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only spherical particles of solid density p s , and use 



vr (Si + s 7 - 



(10) 



and Sj is the assumed radius of a particle of mass rrij. The two are related via m = ^irp s s 3 , where 
p s = lg/cm 3 as in Windmark et al. (2012b). 



Sticking, Bouncing and Fragmenting. The complex dependence of possible collisional out- 
comes on the collisional velocity and particle properties was discussed in Section [TJ and is still very 
much the subject of ongoing investigations. However, since our purpose is not to give quantita- 
tively accurate predictions for the evolution of the particle size distribution function, but rather to 
study qualitatively the impact of the introduction of our new model, we use the simplest possible 



prescription for the fragmentation and sticking probabilities. As in Windmark et al. (2012b) we 
assume that, in any specific collisional event, particles stick for Ajj < vi, and fragment if Ajj > vj. 
Between v b and Vf lies the "bouncing" region. Hence we select, for individual collisions, 

4 = H(Aij - v f ) , 

et j = H(v b -A ij ) , (11) 
where H is a Heaviside function. In all that follows, we use for the sake of comparison with the work 



of Windmark et al. (2012b), the values v b = 5cm/s (unless otherwise specified) and Vf = lOOcm/s. 



Relative velocity. The various possible sources of relative motion for solid particles in proto- 



planetary disks, as reviewed by Weidenschilling & Cuzzi (1993) for instance, can be divided into 



two categories: regular motion caused by frictional interaction with the mean component of the gas 
velocity (usually divided into radial drift, azimuthal drift, and vertical settling), and chaotic mo- 
tion caused by collisions with gas molecules (Brownian motion) and interaction with the turbulent 
component of the gas velocity. 

Let us begin by considering the motion of particles in a gas at rest (i.e. supposing that the 
macroscopic gas velocity is zero in the frame of reference considered). Each particle undergoes 
Brownian motion via collisions with the gas molecules, and thereby acquires a random velocity, 
which is isotropic and has an amplitude drawn from a Maxwellian p.d.f.. Since the velocities of two 
particles both undergoing Brownian motion are uncorrelated, the p.d.f. of their relative velocity 
Ajj is also a Maxwellian with 

f2 A 2 - / A 2 - \ 

m( W;<^(-2<^) ' (12) 

where (vfj) 2 = kT(mi+mj)/mimj, k is the Boltzmann constant and T is the local gas temperature. 
The mean value of the relative collisional velocities is 



xb 8kT(mi + mj) 18 B 

A ^' = V nm imj =yjn a »- (13) 
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Note that, as defined, ofj is exactly 1/V3 times the rms collisional velocity and A.f- is y/8/3ir times 
the rms collisional velocity This implies that for a Maxwellian, to a good approximation, the mean 
and rms velocities are equal to one another - a well-known result that we shall use later. 

In protoplanetary disks, however, particles can have significant mean velocities with respect 
to the gas. This is particularly true for the larger particles, which orbit around the star at 
near-Keplerian speeds, while the azimuthal gas velocity is notably sub-Keplerian (since the gas 
is pressure-supported while the particles are not). As a result of this differential motion, collisions 
with gas molecules transfer net momentum to the particles. This effect is often modeled as an added 
drag term in the particle's equation of motion, and leads to radial and azimuthal drift within the 



disk ( 


Whipple 


1972 


Weidenschilling 


1977, 


Nakagawa et al. 


1986 


), as well as 


toward the mid-plane over time ( 


Goldreich & Ward|1973 


Dubrulle et al.|1995 



Let's consider a frame that is rotating with the local Keplerian angular velocity Ok- In that 
frame, assuming that the surface density of solids is much smaller than the surface density of the 
gas, particles of mass rrtj have mean radial and azimuthal velocities Ui and Vi respectively, with 



1 



1 



(lig - 2SiT]VK) , 



v; 



2S 



(14) 



where u g is the radial gas velocity, and rj is related to the deviation of the gas azimuthal velocity 
from the local azimuthal Keplerian velocity vk = t^k- All these quantities depend on the selected 
disk model, and are presented in more detail in Section |4| Finally, Si is the Stokes number of 
particles of mass riii, defined as 



SiPs 
PmC 



2irsip s 



(15) 



where n is the particle stopping time, p m is the mid-plane disk density, c is the local sound speed 
and S is the local surface density of the gas (see Section [4]). Note that we have assumed here 
that particles remain in the Epstein gas drag regime, which is true for small particles but may 
not be accurate above a certain size. This assumption should be dropped should one require a 
quantitatively more accurate prediction for the evolution of the particle size distribution function, 
but is satisfactory within our stated goals (see Section [T]). 



The mean settling velocity of a particle of mass rrii is calculated as in Birnstiel et al. (2010): 

1 



W: 



-/ijf^K min Si, 



(16) 



where hi is the vertical scaleheight of particles of mass mi, and where the minimum function 
guarantees that the vertical velocity is at most equal to that of the epicyclic motion. The particle 
disk scaleheight hi can be obtained by solving an advection-diffusion equation ( Dubrulle et al.|1995 



-14- 



Garaud||20Q7[ ), which yields 



1 

where h is the local gas disk scaleheight, D{ 



SjflR 1 



(17) 



is the reduced turbulent diffusivity for the 



particles of size Si, and v is the local turbulent viscosity of the gas (see Section [4]). 

As should be obvious from these derivations, particles of different sizes have different velocities 
relative to that of the gas. As a result, they also acquire relative velocities with respect to each 
other. This time, however, no random effects are involved so the mean value of the relative velocities 
of particles i and j, both induced by gas drag, is 



D 



U i 



+ (Vi 



(18) 



and the p.d.f. is well approximated by a 5— function centered around the mean. 



However, the effect of gas drag described above is only strictly valid in a "non-turbulent 
accretion disk" - somewhat of an oxymoron. As first described by Voelk et al. (1980), particles 



interacting with turbulent eddies have inherently stochastic motions. Their velocity p.d.f. is non- 
isotropic (when the turbulence itself is anisotropic), and depends sensitively on the particle stopping 



time compared with the eddy turnover time at the energy dissipation and injection scales (Voelk 



et al. 1980 Ormel & Cuzzi 2007). To add to the complexity of the problem, particle trajectories 
are no longer necessarily independent of one another - two small particles trapped in the same 
eddy have strongly correlated motion. Recent progress has been made to characterize the problem 



using full numerical simulations (Carballido et al. 2010) and simplified vortex gas models (Rast 



& Pinton 2011), although the results still have limited applicability in both cases. A theoretical 



model of the rms relative velocities in turbulence, which we assume is very similar to their mean 
relative velocity Afj by analogy with the case of Brownian motion, was proposed by 



(2007) (see their equations (27)-(29)): 



Ormel & Cuzzi 



A 



l S l - Sj 
Si + Sj 



5 



Si + Re" 1 / 2 Sj + Re" 1 / 2 



1/2 



2 . 2 S{ S j ~\~ 



2S? 



+ 



if Sij < Re" 1/2 

1/2 



3 ' Si + Sj I 2.6 ' IfiSf + SfSj 



if Re" 1/2 < Si < 1 for all j < i , 



1 1 

+ 



l + Si l + Sj 



1/2 



if 1 < Si for all j < i 



(19) 



where Re is the local Reynolds number, which we will take to be 10 8 (see Section B, v e = \J~ac is 
the typical eddy velocity, and a is the turbulent intensity parameter (see Section [4|. Note that in 
these equations we have assumed that Si > Sj. If the opposite is true, then the indices should be 
switched. 
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3. Kernel prescriptions — old vs. new 



All the information described in the previous Section now needs to be combined to create 
the coagulation and fragmentation kernels. In this section, we contrast the "old approach" , which 
includes all ensemble models of particle growth until 2011, and the "new approach" first proposed 
by Galvagni et al. (2011) and Windmark et al. (2012b), which we build upon. 



3.1. Traditional approach 



As discussed in the previous section, the collisional and fragmentation kernels are usually 
given by equation ^ . The problem then reduces to calculating a single value for the mean relative 
velocities of the two particles Ajj, combining information from all the possible sources of motion 
described in Section 



2.2 



Traditionally, A ij is constructed as in Tanaka et al. ( 2005 ) , with 
Ay = JiAf^ + (Ag)* + (A£) 2 , 



I.)' 



and the fragmentation and coagulation probabilities ef ■ and e{- are simply taken to be 



(20) 



' , ^'jJ an d e ij 



(21) 



where the functions ef^ and e\- are constructed (for instance 2 ) as in equation ( jllj . 



3.2. New approach 

In order to take into account the stochastic nature of the particle velocities, collisional and 
fragmentation kernels need to be re-writterj^] as 

/•oo 

Kij = an / Ayp(Aj i )e|-(Ai i )dAy , 
J o 

/>oo 

<i,j J A ii p(A ii )ef j (Ay)dA ij , (22) 

which one may recognize as a much-simplified version of equation ([!]), where p(Aij) is a single 
p.d.f. for the relative velocities of the two particles (dropping the suffix A on the p for clarity of 
notation), which ought to combine information from all the possible sources of motion described 
in Section |2.2| The main difficulty here comes from the fact that, while the p.d.f.s of particles 
undergoing Brownian motion or mean drift are known, we still do not have any knowledge of the 

2 In most prior work to date, more elaborate prescriptions for efj and e{j are used, which include piecewise linear 
or piecewise parabolic functions. 

3 Note that here we have ignored the variability in the area cross-section, as discussed earlier. 
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shape of the p.d.f. of relative velocities for particles undergoing collisions via turbulent motions. 
One is left to chose the latter somewhat arbitrarily. 



Galvagni et al. (2011) and Windmark et al. (2012b) both proposed that A™ should be dis- 



tributed as a Maxwellian. Galvagni et al. (2011 ) considered all possible regimes simultaneously, and 



constructed the Maxwellian so that its mean is given by Ay as defined in equation ( 20 ) . Windmark 



et al. (2012b) only considered the turbulent regime, and required that the mean of the Maxwellian 
be Afj instead. Both models naturally yield the same answer in a regime that is dominated by 
turbulent motions (i.e. when A^- ~ Ajj 3> Af-) but both models fail in the limit where the system 



is dominated by regular motion. The proposed model by Galvagni et al. (2011) overestimates the 



particle dispersion in that limit, while Windmark et al. (2012b) do not address the question. 



In what follows, we build on the work of Okuzumi et al. (2011 ) and propose an alternative, more 



rigorous approach for the construction of p(Aij), which takes into account the fact that individual 
particles have deterministic velocities (induced by radial drift, azimuthal drift and vertical settling), 
in addition to stochastic motions induced by Brownian motion and turbulence. For simplicity, in 
all that follows we assume that the stochastic motions are isotropic]^] 

Let's focus first on a given direction (taking here, for the sake of illustration, the radial direc- 
tion). We model the one-dimensional velocity p.d.f of particle i as a Gaussian with mean Ui (the 
mean radial drift velocity discussed earlier) and standard deviation o~i (specified later). Hence, 



p r (ui) 



2-KO-j 



exp 



Uj 



Ui) 



2a 



(23) 



where the index r denotes the radial direction. Note that, in this context, Ui can be positive or 
negative. Given a second particle indexed by j, with a similarly defined radial velocity p.d.f., a 
well-known (although non-trivial) result from probability theory is that the p.d.f. of their relative 
radial velocities is 



p r \Ui 



Ui 



Pr(A r i 



'27TOV 



exp 



(A 



(Ui - Uj))' 



(24) 



where a?- 



af + a 2 . Note that A f 



sign upon permutation of i and j. 



Ui - 



■Uj can be positive or negative. Furthermore, it changes 



Similar expressions can be obtained for the p.d.f.s of relative velocities in the azimuthal (if) 
and vertical (z) directions. Assuming independence of the distributions in each respective direction, 
the p.d.f. of the three-dimensional relative velocity is simply the product of the three derived p.d.fs: 



p 3D (A ij )dAi 



Pr ( A r> ij )pip (A^jj )p z (A Z ij)dA r ijdA t p^ijdA Z; ij . 



(25) 



4 This assumption is not absolutely necessary, but is very useful. Without it, many of the integrals involved in the 
derivation of the p.d.f of the relative velocities have to be evaluated numerically. 
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From this expression, with a little bit of algebra (see Appendix A) and assuming isotropy in the 
particle stochastic motion (i.e. assuming that cr, and Oj are independent of the direction selected), 
one can show that the p.d.f. of the amplitude of the relative velocity A,,- = 



IS 



2tt<7 - Ag 



(A 



exp 



ij 



A D -) r 

it I 



exp 



(A 



13 



+ Agr 



(26) 



where A® is given by equation (18). This time, both Ajj and Ag are always positive by construc- 



tion, and invariant under permutation of i and j. It can be shown that this p.d.f. is always positive, 
and that the integral over all possible values of A^ is indeed 1, as required. 

In the limit where Ag <C (Tij, which corresponds to cases where the deterministic relative 
velocities of the particles are much smaller than their rms velocities (i.e in the case where relative 
motion is dominated by Brownian motion or turbulent motion), p(Aij) reduces to a Maxwellian 
distribution (to lowest order in a Taylor expansion in Ag), with 



p(A*i) ^ 



2 A 



7T 0~ 



exp 




(27) 



which has a mean 



A 



a 



ij ■ 



(28) 



In order to ensure that the mean velocity of that distribution is the same as that obtained from 
Brownian motion or turbulent motion when these dominate the dynamics of the system, we take 



7T 



O"; 



+ m 2 



(29) 



1-3 



where Ag and Af^ are defined in equations (13) and (19) respectively. Our formalism recovers the 



idea proposed by Galvagni et al. (2011) and by Windmark et al. (2012b) in the limit where turbu- 



lence is the dominant factor. Note that, as defined, erjj is not the rms velocity of the distribution 



defined in equation (26), although it is related to its variance. In the limit where the mean drift 



velocity is close to zero, <Jij is, as discussed earlier, l/\/3 times the rms velocity. Our derivation is 



so far very similar to that of Okuzumi et al. (2011), and recovers their relative velocity p.d.f. in 



the limit where the particle dispersion is dominated by Brownian motion. 

When Ag S> (Tij, which correspond to cases where the deterministic relative velocities of the 



particles are much larger than their rms velocities, the second Gaussian term in equation (26) is 
negligible, and the relative velocity p.d.f. reduces to 



P(Ay) ^ 



1 A 



27TCTy A^- 



P exp 



(Ay 



Ag) 2 ' 



(30) 



It it interesting to note that this is not a Gaussian, but has a somewhat more significant tail for 
larger values of A^. 
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It is worth noting that the mean relative drift velocity, A?, is not equal to the actual mean 



relative velocity Ay; the latter is given by 



A 



ij 



Ayp(Ay)dA 



A D - 



' A D - 

\/2ai. 



+ 



cry exp 



(A 



D\2 
ij) 



(31) 



The relationship between Af- and Ay is shown in Figure [lj We see that Ay tends to Af- when 
A® 3> aij, i.e. when the mean particle velocity is large compared with the rms particle velocity, 
and to \J 8 /ira^ when Af- <C cry, i.e. when the mean particle velocity is small compared with the 
rms particle velocity. This behavior is consistent with expectations. 




Fig. 1. — Comparison between the mean particle velocity and the mean drift velocity. The solid line 
shows Ay /aij. The slanted dotted line is the y = x line, showing that Ay — > Af- when Afj 2> cry. 
The horizontal line is at Ay = y/8/ ira^ . 



Finally, we can now calculate the new coagulation and fragmentation kernels, following equa- 



tion (22), by convolving the velocity p.d.f. constructed in equation (26) with the individual collision 



sticking and fragmentation probabilities. With the simple expressions for e\- and e\- given by equa- 
tion (11), it is possible to calculate these kernels analyticalljj^] 



Fij = aij J Ayp(Ay)e^-(Ay)dAy = ay J Ayp(Ay)dA 



f (A^) 2 + cr 2 . 

- n -AlVt) 



a; 



1 



2-A^ 



and 



poo rv b 

/ Ayp(Ay)ey (Ay )dAy = ay / Ayp(Ay)dAj 
Jo Jo 



(32) 



(33) 



5 For other more complex prescriptions for e*'f , the integrals need to be performed numerically. 
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ij 



'■>.) 



2A^ 



2erf 



' A fl 



A(« 6 ) 



+ 



2A£ exp 




where for simplicity of notation we have defined 



( v + A D ^ 
A{v) = erf ' 



V2a 



( v - A^ v 
erf 1 



>/2. 



= (u + Ag) exp 



(v-A°f\ f (v + A?)* 



(v - A£) exp 



2 4 



Note that if A^ <C then (to lowest order in a Taylor expansion in A5) 



2 4 + j 
2 4 



exp 



A' 



'J 



Aij 



1 24+^ 
1-^2— exp 




(34) 



(35) 



If desired, one can thereby construct "mean" sticking and fragmentation probabilities 



Kij Xjl>(Xj)<1j(Xj)<IX.i 



J °° Aijp(Aij)dAij 
AijpiA^A^dA 



/ °° A^A^dA 



(36) 



so that, using these expressions, we have an exact analog for equation ([9]). Note that the second 
term in each of these expressions is a rather general formula for the mean sticking and fragmentation 
probabilities, which is valid for any assumed velocity p.d.f., and any assumed individual collision 
and fragmentation probabilities, but does assume that the collisional outcomes are independent of 
the particle density and shape. 



It can easily be verified that when efj(Ajj) + e{-{Aij) 



1 for individual particles, then we also 



have efj + e£ 



1, or equivalently, 



K-ij -\- Fij — QijAij . (^^) 

This happens for instance when Vb = Vf (i.e. in the absence of bouncing) in the simple Heaviside 



model given by equation (11), and can be verified directly by adding equations (32) and (33). 



Note that the coagulation kernel derived in equation (33) looks like the one presented by 



Okuzumi et al. (2011 ) in their equation (12), but is slightly different because the sticking probability 
of individual collisions in their paper is taken to be a piecewise linear function of the collisional 
energy rather than a Heaviside function of the collisional velocity. Our results would of course 
recover one another had we chosen the same sticking probability function, and in the limit where 
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the particle dispersion is not affected by turbulence (i.e. in the limit where Af 



0). By contrast, 

Windmark et al. (2012b) would not recover the same formula in any limit. Indeed, note that 



Windmark et al.| (2012b) discuss the role of the total sticking and fragmentation velocities (see 
their equations (4) and (5)), in our notation written as 



i e ij] = / PM(Aij)dAij and [e{] = I p M (^ij)dAij , 

Jo Jv s 



(38) 



(where pm denotes a Maxwellian distribution) rather than the mean sticking and fragmentation 
velocities efj and e{-. However, as should be clear from the derivation presented above, these "total" 
probabilities play no role in the calculation of the kernels - the means must be used instead. It is 



not clear from their discussion whether Windmark et al. (2012b) used these expressions directly to 



calculate the kernels (since the latter are not explicitly written out), but if they were, then their 
numerical results should be considered incorrect even though they look qualitatively similar to ours 
(see Section [5]). 



Figure [i] shows ef^ and e\- , as calculated in equation ( 36 ) , as a function of cr^ 



>'b 



and Afj, for 



5cm/s and vj = lOOcm/s (see Section 2.2). We see that when both o~ij and Af- are small 



compared with the bouncing threshold, sticking is very probable, while when either o~ij or Ajj are 
large compared with Vf, fragmentation is very probable. The transition between the two regimes, 
however, is much smoother when o~ij is large than when it is small; this behavior is, again, consistent 
with expectations. 

It is important to note that the effective sticking and fragmentation probabilities efj and 
are never zero, by contrast with the individual collision functions efj and e{j. This result is very 
general, and stems from the fact that eiy and e{- are convolution integrals of e?- and e{- with the 
assumed particle velocity p.d.f.s times the particle relative velocity. As a result, fragmentation 
is always possible (albeit unlikely) even when the mean relative velocity of two particles is very 
small, and sticking is always possible (albeit unlikely) even when the mean relative velocity of two 
particles is very large. This has a number of important consequences, as shown in Section [5} 



4. Disk model and consequences for particle dynamics 

Since the mean and stochastic relative velocities of the particles depend on the disk model and 
the location considered we present them here for completeness, and discuss them in the light of the 
previous section. 
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Fig. 2. — Effective sticking and fragmentation probabilities, assuming a bouncing threshold of 
Vb = 5cm/s and a fragmentation threshold of vt = lOOcm/s, and calculated using equations (36). 
The left-hand plot shows shows efj and the right-hand plot shows qj-. The color bar is the same 
for both plots. Note how perfect sticking requires that both a^j and should be smaller than uj,. 



Also note how the transition from perfect sticking to perfect fragmentation, across the "bouncing 
region", is much smoother when increasing cr^ at fixed A® than increasing Kf- at fixed aij. The 



Maxwellian-only approach (cf. Galvagni et al. 2011 and Windmark et al. 2012b) assumes that 
A D - = 0. 
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4.1. Disk model 



We assume that the surface density of the gas follows a truncated power law: 

M 



-r/R 



(39) 



2irRr 

where M is the disk mass, r is the orbital radius, and R is the truncation radius. As discussed 



by Lynden-Bell & Pringle (1974) (see also Hartmann et al. 1998 Garaud 2007), this model has 



the desirable property of being an attracting self-similar solution of the equations describing the 
"viscous spreading" of the disk, as long as the turbulent viscosity v varies linearly with radius: 

r 



u{r) = I/AU: 



(40) 



1AU ' 

where the subscript AU refers to a quantity measured at 1 AU. Conservation of angular momentum 
implies that the radial velocity of the gas is 



u g (r) 







(r 1 / 2 i/(r)S(r)) 



(41) 



rV2E(r) dr V " v ' v ') "1AU V2 
while the radial force balance near the mid-plane implies that the azimuthal velocity of the gas is 

1 dp r , 



v g {r) = v K (r) + - / s n / n o 
2p m (r)tt K (r) dr 

where p m is the mid-plane gas pressure, and where 

1 

7 ?( r ) = 



(1 - rj(r))v K (r) 



2rp m (r)ti^{r) dr 



(42) 



(43) 



As illustrated by equations (41) and (42), the local gas velocity depends on the local thermody- 



namical structure of the disk. We now turn to describing the latter in more detail. 

Assuming that the disk is thin, vertically isothermal and in hydrostatic equilibrium implies 
that the density profile is 



p(r,z) = p m (r)exp 



where the disk scaleheight h(r) is related to the local sound speed c(r) via 

h{r ) = 

Integrating the density profile across the disk yields 

Pm(r) ■ 



(44) 



(45) 



(46) 



Meanwhile, the mid-plane gas pressure is given by the equation of state (which we assume here to 
be a perfect gas), so 

7? n..Jr\T(r\ 

(47) 



= n Pm (r)T(r) = c2{r)pm{r) 



- 23 - 



where is the mean molecular weight of the gas and 1Z is the gas constant. 

In both cases, we need to know the local sound speed to evaluate p m and p m . Following the 
standard a— model for turbulent viscosity, we use 

v{r) = ac(r)h(r) , (48) 



where a is usually assumed to be constant in the entire disk. Combining equations (40) and (|48 
we have 



2( \ 2 / r \~ 1/2 i 2 ^AU, 



c (r) = c AU [jj^) where c AU = — Ouj , 

iau J where ^ AU = ' (49) 

Having laid out the governing equations for the disk, we need to choose its actual parameters. 



For the sake of comparison with the work of Windmark et al. (2012b), we select the parameters for 
our disk to have the same properties at IAU (see their Table 1). This immediately yields a = 10 -4 , 
cau = 10 5 cm/s, and the Reynolds number Re = 10 8 . Using a mean molecular weight of [i = 2.3, 



the temperature at IAU is 276K, close to the value reported by Windmark et al. (2012b). They 
choose a gas surface density of £ = 1700g/cm 2 , which we can recover to a good approximation 
by selecting (for instanc^]), a central star of mass M+ = 0.75Mq, surrounded by a disk of mass 
Md = 0.05M*, with a R = 30AU truncation radius. Using these parameters we find that the viscous 
timescale at IAU is about 400,000 years, so the gas disk would not evolve significantly during the 
simulation. 

At 30AU, our disk model has a local gas surface density of 21g/cm 2 and a temperature of 50K. 
At this radius, the viscous evolution timescale is larger than lOMyr, so again, the gas disk would 
not evolve significantly during the simulation. Finally, in all that follows we select a dust-to-gas 
mass ratio of Z = 0.01. This ensures that the dust dynamics cannot influence the gas dynamics, 
as assumed throughout this work. 



4.2. Particle dynamics in the disk considered. 



In what follows, we shall be interested in two representative disk regions, located close and 
far from the central star respectively. The close-in region is selected to be at IAU, for ease of 
comparison with the work of Windmark et al. (2012b), and the far-out region is selected to be at 
30AU, to address the problem of particle growth beyond cm-size discussed in Section [T] 



6 Unfortunately, Windmark et al. (2012b I do not report on the mass of the star used in their simulations; since the 
stellar mass determines the dynamical timescale of the disk, precise comparisons between our simulations and theirs 
is impossible. 



-24- 



As shown in Figure [3j relative particle motions are dominated by turbulence for sizes up to 
about 10 cm at 1AU. Collisions induced by differential radial drift dominate for larger particles. 



Hence, we expect our model to become more relevant than the one proposed by Galvagni et al. 



( |201l[ ) and |Windmark et al.| ( |2012b[ ) if growth occurs beyond 10cm at lAUQ At 30AU, this tran- 
sition happens for mm to cm size particles; the model proposed here should therefore be used 
preferentially to model growth beyond mm-size, which are indeed the ones we are interested in. 

Figure [3] also shows the logarithm of the mean relative velocity Ay of all particle pairs, 
and contours of the bouncing and fragmentation thresholds Vb = 5cm/s and Vf = lOOcm/s re- 
spectively. In what follows, we shall often refer to these contours as the "original" bouncing and 
fragmentation barriers, as they often play that role in simulations which use the traditional method 
for kernel construction (see Section 3.1). In any case, we find that bouncing is expected to affect 
particle growth beyond mm-size at 1AU, and 10 micron-size at 30AU, for the selected value of V&. 
Fragmentation is expected to be important for particles around cm-size at 1AU, and mm-size at 
30AU. 



5. Results 

We now proceed to study the collisional evolution of the particle distribution function, inte- 
grating equations ^ numerically at 1AU and 30AU in the disk, and compare the results obtained 
using three possible kernel construction procedures: 



Model O: The "Old" approach, in which kernels are constructed using equations (|9|), with 



Aj,- given by equation (|20l) and ef - and q - given by equation (|21). 



Model M: The "Maxwellian" approach, in which kernels are constructed using equation (35), 
with Ajj given by equation (28) and cry given by equation (29). This recovers in spirit 



(although not in detail) the idea proposed by Galvagni et al. 



2011), and Windmark et al. 



(2012b). 



• Model N: The "New" approach, in which kernels are constructed using equations (32) and 



(33). This treats mean motions and stochastic motions separately, as discussed in Section 3.2 



In all cases, the mean area cross section is given by equation (10). 



In all simulations shown, the initial conditions are nearly mono-disperse, and are taken to 
be a narrow Gaussian centered on a mass-point that is slightly larger than the minimum mass- 
point considered. They have negligible influence on the results, except in the case of Model O 



7 Note that as radial drift becomes more important, the local model approximation unfortunately begins to fail; 
in this sense, the results of our model should only be considered indicative of processes which should be taken into 
account, rather than actual size predictions. See Section IB] for detail. 
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Fig. 3. — Collisional regimes and collisional velocities (top figure, 1AU, lower figure, 30AU). In 
each plot the top triangle shows the dominant collisional regime (purple = turbulence, orange = 
radial velocity, yellow = azimuthal velocity). Note that, for the disk regions and particle sizes 
considered, Brownian motion and vertical settling never appear to dominate. The lower triangle 



shows the logarithm of the actual mean relative velocity Ay as calculated in equation (31), with 



the corresponding color bar shown on the right, as well as contours representing the bouncing and 
fragmentation thresholds: Ay = 5cm/s (left) and Ay = lOOcm/s (right) respectively. 
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when bouncing is included (see below for detail). Unless otherwise specified, the mass- mesh is 
logarithmically distributed and contains 6 points per decade in mass. In what follows we plot the 
surface density distribution function 



evolves towards a power-law with 
proportional to m 2 ~\ 



dn 
dm 



dm 



hm 2 ^. As a result, if the mass distribution function 
oc m~ x , then the surface density distribution function is 



At 1AU, as in Windmark et al. (2012b), we begin by ignoring the possibility of bouncing, 



then include it. We finish by considering a model in which high-mass-ratio collisions lead to mass 
transfer from the projectile to the target ( Teiser <fe Wurm|2009 Kothe et al.|2010 ), as in Windmark 



et al. (2012a) and Windmark et al. (2012b). We then use the latter model to investigate the time 



evolution of the particle distribution function at 30AU as well. 



5.1. No bouncing, 1AU 

Let us first look at particle growth at 1AU, in the case where the possibility of bouncing is 
ignored (taking Vb = vt = lOOcm/s in all simulations). We compare the three possible evolution 
scenarios laid out above. We find that in all cases the particle mass distribution function rapidly 
evolves to a steady state (within less than 5000 yrs), shown in Figure |4j In each case, the steady 
state distribution dn/dm is well-approximated by a truncated power law, with the same slope as 



that of the fragment mass distribution (see Section 2.1). The truncation size, however, depends on 
the model considered. 

In Model O, particles grow up to a few centimeters in size. Predicting the maximum mass/size 
achievable in this model is in fact very easy - one can simply read it at the intersection of the 
horizontal axis and the Aij = vj contou^ in Figure [3j Indeed, since smaller particles are much 
more numerous than larger ones, high- mass-ratio collisions are much more common than equal- 
mass collisions. The growth of a large particle "i" is thus stalled when collisions with much smaller 
ones begin to lead to fragmentation instead of sticking. In Model O, where this transition occurs 
through a Heaviside function, the maximum particle size Si achievable is then simply given by 
Ajj = Vf, taking j = 1 to identify collisions with the smallest particles in the simulation. 



As expected from the discussion in Section 3.2, since particles are not seen to grow up to a 



size where radial drift becomes important, the "Maxwellian" (M) model and the "New" (N) model 
predict very similar outcomes. In both cases, the maximum size achieved is significantly smaller 



than in Model O. This is quite different from the results of Windmark et al. (2012b), and could 
be due to their possible use of "total" probabilities instead of the mean fragmentation and sticking 
probabilities (see Sections. Furthermore, we also find that the maximum particle size achievable 



technically, since this method only concerns Model O, the maximum size should be read from an equivalent plot 



showing the contour of Ajj = i>/ as calculated through equation (201 rather than (311; however, it happens that the 
two contours look very similar in the region considered. 
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is sensitive to the minimum particle size considered in the simulation. This is illustrated in Figure 
[4] which shows, in addition to the three models described above, the outcome of a different run 
using Model N, with exactly the same parameters and same resolution, but simply decreasing the 
minimum size considered by two orders of magnitude. The maximum particle size achieved in this 
case is three times smaller than before. 

After further investigation, we found that this effect is in fact generic to any model for which the 
fragmentation probability does not drop exactly to zero below a certain threshold velocity (which 
is the case in Models M and N, but not in Model O) - and is in fact very easy to understand from 
a physical point of view. Indeed, as discussed above, collisions are much more frequent the smaller 
the projectile. In the very crude collisional model considered in this Section, it is possible (albeit 
unlikely) for a micron-size particle to collide with a much larger one and cause its fragmentation. 
Even though the fragmentation probability is very low, the sheer number of collisions that take 
place ensure that fragmentation does happen on a regular basis. Since the number density of 
particles is a steep function of their mass (or size) , the new maximum particle size achievable thus 
sensitively depends on the smallest particle considered in the numerical integration, as well as the 
rate at which the fragmentation probability drops to zero when Ay drops below Vf. This can be 
shown analytically, and will be the subject of a separate publication. 

We note, however, that this trend disappears when considering models in which high-mass-ratio 



collisions do not lead to complete fragmentation, see Section 5.3 for details 



5.2. With Bouncing, 1AU 

We now consider the effect of bouncing, by setting Vb = 5cm/s. Again, we compare the three 
scenarios discussed above. We set the minimum particle size in the simulation to be Smin — 0.01cm, 



for ease of comparison with Windmark et al. (2012b) and with the results of the previous section 



bearing in mind the potential effect of that choice on the outcome. In all cases, the simulation is 
integrated for 30,000 yrs, and the resulting particle size/mass distribution functions are shown in 
Figure [5} 

At that point in time, Model O is still evolving. By construction, it will continue to do so 
until all particles have grown to a size above which any collisions they may undergo only result in 
bouncing. Little by little, all small particles disappear from the system. In many ways, this is the 
worst-case-scenario in terms of comparison with observations: small particles are entirely depleted, 
but growth to sizes larger than the bouncing barrier is also impossible. Note that the outcome 
in this case is strongly dependent on the initial conditions chosen. Suppose for instance that all 
particles are initialized within the bouncing region - then, the particle size distribution function is 
in fact in a steady state. Here, we have selected to initialize all particles within the sticking region, 
to emphasize the effect of the bouncing barrier. 



By contrast with Model O, Models M and N have already reached a steady state by 30,000yrs. 
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Fig. 4. — Particle evolution at 1AU, after 5,000 years, in models where the bouncing and fragmen- 
tation velocities are Vj, = vj = lOOcm/s. Models O, N and M are described in the main text. In 
these simulations, the size- mesh spans the range s S [10 -2 , 10 4 ]cm (although the figure is truncated 
to focus on the region of interest). Model N (2) is the same as model N, but with a numerical mesh 
that extends down to a minimum particle size of 10 _4 cm instead of 10~ 2 cm, with the upper mass 
point correspondingly reduced to 10 2 cm to keep the same resolution. 



As in the previous section, the two are very similar to one another, as expected from the fact 
that particles do not appear to grow beyond a size for which radial drift becomes important. As 



discussed by Windmark et al. (2012b), the maximum particle size achieved is, this time, larger than 
in Model O since occasional low-velocity collisions are still possible even when the mean relative 
velocity of the two particles is already beyond the bouncing barrier. More importantly, although 
fragmentation is rare (since the mean fragmentation probability, at the sizes achieved, is very small), 
it is nevertheless sufficient to replenish the small-particle population. In other words, models M 
and N provide more growth than Model O, and maintain a significant population of small grains, 
but the largest particle size achieved at steady-state is still much smaller than what is needed to 
form planetesimals. 



5.3. With Bouncing and mass transfer, 1AU 

As discussed in Section [TJ high-mass-ratio collisions are unlikely to result in the complete 



destruction of the larger body, as assumed in Sections 5.1 and |5.2| Instead, the smaller particle 
(dubbed the "projectile" hereafter) is more likely to excavate a small crater from the larger one 
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Fig. 5. — Particle evolution at 1AU, after 30,000 years, in models where the bouncing velocity 
is reduced to Vf, = 5cm/s compared with Figure |4j The size range considered in Model O was 
reduced to s E [10 -2 , l]cm, but with increased resolution (with 20 mass-points per decade in mass) 
to capture the accumulation of particles near the bouncing barrier. Models M and N, by contrast 
have the same resolution and size range as in Figure [4] (with s E [10~ 2 , 10 4 ]cm and 6 mass points 
per decade in mass), even though we only show sizes up to 10cm in this plot. 



(dubbed "the target" hereafter) and/or could itself partially or completely stick to the target 
( |Langkowski et al.|[2008| |Teiser fc Wurm||2009| |Kothe et aLl|2010[ ) thus resulting in mass transfer. 

Unfortunately, our understanding of the conditions under which mass transfer occurs, and 
in particular its dependence on the particle masses, porosity and velocities, is still in its infancy. 
For this reason, we use here a very simplified model of the process, in which we assume that the 
projectile completely sticks to the target if the mass ratio of two particles is larger than a certain 
critical value (ft (which we take to be 50, by analogy with Windmark et al. (2012b)), and if the 
collision would have otherwise led to fragmentation. To implement this numerically, one simply 
needs to use the following algorithm: (1) Calculate the kernels K^ and F^ in the absence of 

first redefine Kij as 



"mass transfer". (2) If the mass ratio of the two particles is greater than 



K 



K 



y + Fij , and then set Fij 



0. 



This prescription differs from that of Windmark et al. (2012b), who assume that only 10% 

% fragments. We recognize 



of the mass of the projectile sticks to the target, and the other 
that our own choice is likely to overestimate the actual amount of mass transfer occurring during 
a collision, but since the actual value of (ft is very poorly constrained anyway, the error made is 
within the same order of the approximation. We also show in Section 5.5 that the results presented 
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below still hold qualitatively with lower mass-transfer rates. The advantage of this approach is 
that the mathematical structure of the effective fragmentation/sticking probabilities is now very 
easy to visualize. For low mass ratio, they are as depicted in Figure |2j For high mass ratio, e\- = 
while efj now looks like the sum of the previously calculated sticking and fragmentation probability 
functions (i.e. the sum of the two surfaces shown in Figure [2j). In other words, efj is close to unity 
for low collisional velocities and for high collisional velocities, with a substantial "gap" (or perhaps 
a "moat" may be more appropriate) in-between that corresponds to the bouncing region. The 



latter is delimited by the original bouncing and fragmentation barriers (see Section 4.2). 



Mass transfer has fundamental implications for particle growth: as long as e|„- never strictly 



drops to zero anywhere, which is the case in Models M and N as discussed in Section 3.2 there 
is always a possibility for growth across the "moat" in high-mass-ratio collisions. Indeed, While 
bouncing in this regime may be the most likely outcome of a collision, there is always a possibility 
of very low- and very high-velocity events that result in coagulation and growth. Hence, once 
particles are large enough, they may always continue to grow by sweeping up smaller ones. As long 
as fragmentation via low-mass-ratio collisions is rare enough, the net effect is growth beyond the 
original barriers. 

Figure [6] presents the results of simulations that include the simple mass transfer model dis- 
cussed above. We compare simulations using Models O, M and N, this time after 10,000 yrs only. 
The effect we are interested in is already very clear, and only becomes more pronounced beyond 
that time. One immediately notices that mass transfer has no effect on Model O. This was already 



pointed out by Windmark et al. (2012b), and is an obvious result in the light of the discussion 



above. Indeed, in Model O, both efj and e{j strictly drop to zero when Ay G [vb,Vf]. This means 
that, even for high mass-ratio collisions, efj strictly drops to zero for Ay > Vb, preventing any 
possibility for growth beyond the original bouncing barrier. 

For Models M and N, however, we see a dramatic increase in particle growth. Furthermore, 
since the largest particles are now in a regime that is dominated by radial drift rather than turbu- 
lence, the predicted evolution of the size distribution functions for models M and N differ signifi- 
cantly. After only 10,000 years, the maximum particle size has increased by more than 4 orders of 
magnitude for Model M (corresponding to more than 12 orders of magnitude increase in mass, as 



found by Windmark et al. (2012b)). The effect is even more pronounced for Model N, where the 
maximum particle size is yet another order of magnitude larger. Furthermore, the surface density 
of larger particles is 3-4 orders of magnitude larger in model N than in Model M. In short, correctly 
modeling the difference between regular and stochastic processes when constructing the particle 
relative velocity distribution function turns out to be highly beneficial to growth. Finally, note that 
even though fragmentation no longer completely suppresses the growth of large particles, it is still 
amply sufficient to replenish the small grain population. 
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Fig. 6. — Particle evolution at 1AU after 10,000 yrs, with bouncing and mass transfer, in models 
where the bouncing velocity is Vf, = 5cm/s, and where the mass ratio above which every collision 
(that would otherwise lead to fragmentation) sticks is = 50. The size range and resolution 
for Model O is the same as in Figure [5} The size range for Models M and N was increased to 
s G [10 -2 , 10 7 ]cm, with a corresponding increase in the number of mass-points to keep the same 
resolution. 



5.4. Evolution of the particle distribution function at 30 AU 

Using Model N, we can now study other regions of the disk. Here, we consider the disk at 
30 AU, and include both bouncing and mass transfer with the same parameters as in the previous 
sections: = 5cm/s, vj = lOOcm/s, and <fi = 50. As seen in Figure[3j the bouncing threshold now 
occurs for much smaller particles than at 1AU. For this reason, we shift the mass-mesh so that 
s £ [10 -4 , 10 5 ]cm instead, and keep the same resolution. 

In order to present complementary and comparable information to that of the previous sections, 
we show the actual temporal evolution of the particle distribution function. It is clear from Figure [7] 
that the system now contains two interacting but distinct particle populations. The small particle 
population is distributed as a truncated power law, with a power index dictated by the assumed 



fragmentation law of the model (here, with dn/dm ~ m~ L83 , see Section 2.1 ), and a maximum size 
around 0.1 mm. At early times, the bouncing barrier induces a pile-up just above that size, giving 
the appearance of a fairly mono-disperse population. However, the pile-up gradually disappears 
later on, as the fragmentation of an increasing number of larger particles replenishes the small 
particle population. Meanwhile, the larger particles steadily sweep up the smaller ones, resulting 
in the gradual emergence of a large-particle population. 
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Fig. 7. — Particle evolution at 30AU with bouncing and mass transfer using model N, where 
the bouncing velocity is Vb = 5cm/s, and where the mass ratio above which every collision (that 
would otherwise lead to fragmentation) sticks is <j> = 50. This figure shows snapshot of the mass 
distribution function at different times. Note the gradual emergence of a large particle population, 
of sizes up to a few cm. Although difficult to see, mass is indeed conserved in this simulation - as 
time evolves, the particle peak is slightly eroded. 
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After about 500,000 years, the large particle population has become very significant indeed. 
Its distribution, especially at later times, can also be viewed as a truncated power law. The power 
index is shallower than that of the small particles, with dn/dm ~ m -L2 . The maximum particle 
size within the large-particle population initially grows quickly with time, increasing by 3 orders 
of magnitude in the first 100,000 years. Later on, the growth of the maximum particle size slows 
down in favor of a steady increase in the number density of the large particles while keeping the 
same overall shape of the distribution. 

It is clear from this simulation that the effects described here have the potential of solving 
simultaneously the two important puzzles raised by observations of protoplanetary disks - the per- 
sistence of a small-particle population for millions of years, and the inferred presence of rather large 
particles (>cm size) far out in the disk. Furthermore, the emergence of two particle populations 
with different mass distribution functions, rather than a single continuous power law, is consistent 



with the observations of Wilner et al. (2005). We will return to these points in Section 6.2 At this 



point, however, it is time to look into the dependence of the model results on the input parameters. 



5.5. Dependence on parameters 

In all previous simulations, we used the same values for the fragmentation threshold Vf, the 
bouncing threshold Vb, and the mass ratio above which mass transfer happens. We now vary the 
latter two to see their effects on evolution of the particle distribution function. The results are 
presented in Figure [8j 

We consider here the same disk at 30 AU (very similar results apply at 1AU). In all simulations, 
we fix Vf = lOOcm/s, but vary Vb between 5cm/s and 20cm/s, and <fi between 50 and 200. We evolve 
the system of equations ^ for 500,000 years. In all cases we observe the same qualitative evolution 
of the particle size distribution function into the two particle populations described in the previous 
section (small particles and large particles). However, we also see that the transfer of mass from the 
small to the large particle population, which in turn controls the large particle population growth 
rate, depends quite sensitively on cfi, and on the width of the bouncing region (i.e. on vj — v\,). 

First, we find that the larger 4> is, the slower the mass transfer from small to large particles. 
This effect is intuitive: a given target mass m, will effectively sweep all particles of mass mi /(f) or 
smaller. The total amount of mass available for growth can be evaluated from $™ % m(dn/ dm)dm, 
and is naturally smaller when cf) increases (the details of exactly how much smaller depends on the 
specific shape of the mass distribution function for small particles). 

The mass transfer rate is also dependent on the width of the bouncing gap. To understand this, 
let's consider for instance the tiniest particles in the system (taking i = 1), and look at their sticking 
probability with other particles around. In a mass transfer scenario, collisions with anything of 
mass rrij > 4>m\ do not lead to fragmentation. For (f> = 50 and s% = 10 _4 cm, this corresponds to 
Sj > 4 x 10 _4 cm - in other words, very few of the collisions involving si lead to fragmentation. 
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However, the mass transfer model considered here only leads to sticking for collisions that would 
otherwise lead to fragmentation. Bouncing remains a barrier to growth, and thus the larger the 
bouncing region, the stronger this barrier is. Figure [9] illustrates this clearly, by showing the mean 
sticking probability of a particle of index i = 1 colliding with a particle of size Sj. The wider 
the gap, the deeper the minimum in the function efj. This minimum is the growth bottleneck of 
the system, and thus controls the transfer of mass from the small particle population to the large 
particle population. 



6. Discussion and conclusions 

In the previous sections, we presented a new model for the construction of the coagulation 
and fragmentation kernels involved in the evolution of the particle size distribution function, which 
takes into account the p.d.f.s of the relative velocities of the particles. In particular, we include a 
more physically motivated approach which separates the deterministic from the stochastic relative 
particle velocities. When combined with the possibility that high-mass-ratio collisions can lead to 
the growth of the (high-mass) target rather than its partial or complete destruction, this model is 
potentially able to solve three of the longest-standing problems raised by observations of protoplan- 
etary disks: the persistence of small-grains for millions of years despite ongoing coagulation, the 
existence of cm-size particles far out in the disk, and finally, the formation of large planetesimals 
close to the central star inferred from the ubiquity of extrasolar planets. 

In this section, we now discuss the model in more detail, and in particular, which of these results 
are model-dependent, and which are not. We then review our result in the light of observations. 



6.1. Model dependence 

The results of Section [5] suggest that growth beyond the traditional "bouncing" and "frag- 
mentation" barriers, and the maintenance of a small-particle population, merely requires two in- 
gredients: (1) a threshold above which high-mass-ratio collisions lead to sticking rather than to 
fragmentation, and (2) mean sticking and fragmentation probabilities that never strictly drop to 
zero. 

In fact, these ingredients are always naturally expected to arise in any "physically motivated" 
model for the evolution of the particle size distribution function. Their absence from previous work 
can only be viewed as unfortunate and misleading oversimplifications of the problem. 



Let us begin by considering the case of high-mass-ratio collisions. In the original work of Dulle 



mond & Dominik (2005) for instance, the fragmentation and sticking probabilities were constructed 
by considering the total kinetic energy of the collision in comparison with the binding energy of 
the particle. In these kinds of models, high-mass-ratio collisions are much less likely to lead to 
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Fig. 8. — Particle evolution at 30AU with bouncing and mass transfer, in models with varying 
bouncing velocity and mass ratio above which every collision sticks, <p. This figure shows snapshot 
of the mass distribution function at 500,000yrs. Note that all models are still evolving at the given 
time. 
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Fig. 9. — Effective sticking probability ef ■ of a particle of index i = 1 (the smallest possible particle) 
with a particle of size Sj, for <fi = 50, Vf = lOOcm/s and various values of Vb at 30AU. The wider 
the bouncing region, the deeper the minimum in the functio n ef ? -. The apparent discontinuity for 
sj ~ 10~ 3 cm is due to the discontinuity in the Ormel &: Cuzzi (2007) turbulent velocity prescription. 
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the fragmentation of the larger body, and sticking becomes a more plausible outcome. In some 
sense, such prescriptions already incorporate the idea of mass transfer, without any need to define 
it artificially. Unfortunately, much of the work published between 2005 and 2011 dropped this 
energy-dependent view of the fragmentation/sticking thresholds in favor of a velocity-dependent 
one (such as the one used, for the sake of illustration, in this paper). We advocate that future work 
should return to using the energy-based approach in which the notion of mass transfer naturally 
emerges, in conjunction with the use of velocity p.d.f.s for the construction of the kernels. 

Mass transfer on its own, however, is not sufficient: the second condition is equally important. 



This was illustrated in Section 5.3, where we showed that even with mass transfer Model O does not 

lead to growth beyond the fragmentation barrier. In other words, in any model where there exists 

a given particle size Sj for which the mean sticking efficiency ef • = for all possible particle pairs 

growth beyond Si is naturally impossible. This was a rather common outcome of models that 

— — f 

used the "traditional approach" in conjunction with piecewise defined functions efj and q- that 



drop strictly to zero across a particular velocity threshold, as in the work of Brauer et al. (2008) 



and the piecewise linear model used by Birnstiel et al. (2010). 



In reality, however, the mean sticking probability for a given particle size is never expected 
to drop to zero entirely. Even when the mean collisional velocity/energy of a particle pair 
is high, low velocity /energy collisions are always possible. This is explicitly taken into account 
when using p.d.f.s of relative velocities in the calculation of the mean sticking and fragmentation 



probabilities (see Section 3.2). Indeed, since the latter can be written as convolution products of 
the relative velocity p.d.f. and the sticking and fragmentation probabilities of individual collisions, 
they are always strictly positive. Again, we advocate that such a model should always be used in 
the future. 

As long as it is using the two ingredients listed above, we expect that any local model will 
yield answers that are similar to the ones shown in Section [5] for Model N, and will reveal the 
emergence of two populations of particles: a small-particle population constantly replenished by 
the fragmentation of larger bodies (with a size distribution function controlled by a collisional 
fragmentation cascade), and a large-particle population that grows by sweeping the smaller particles 
(with a size distribution controlled by a coagulation/fragmentation balance). In this sense, anyone's 
preferred prescription for the relative velocity p.d.f.s, for the mass transfer scenario, for the effect 
of porosity on the collisions, should yield qualitatively similar results - the latter are certainly not 
model-dependent. 

However, we also showed that the details of the mass transfer rate between the two populations 
are, unfortunately, very much dependent on the model considered, and within a given choice of 



model, on the selected parameters. We illustrated this in Section 5.5 by considering a range of 
parameters within our own toy model, and showed that the predicted surface density of large 
particles in the system can vary by 5-10 orders of magnitude simply by changing the parameters by 
a few. This somewhat unsatisfactory result can be understood, as shown earlier, by the fact that 
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the mass exchange rate between the small and large particle populations is ultimately controlled 
by the sticking efficiency bottleneck, as well as the available surface density of small particles that 
a larger one can "sweep", which are themselves strongly dependent on the parameters. 

As discussed in Section [TJ the model presented in this work was, by and large, selected for 
simplicity rather than physical realism. Comparing results directly with observations will require 
much more sophistication in the physics included. The particle structures (porosities, composition) 
should be taken into account. The sticking and fragmentation probabilities for a given collision are 
more likely to depend, in a complex manner, on the collisional energy and the material strength 
of the two particles. This should also be taken into account instead of using equation (11). In 



addition, the relative velocity p.d.f.s of two particles undergoing turbulent motion is not likely 
to be a Maxwellian, as we had to assume here. Indeed, the velocity p.d.f. of a single particle 
interacting with turbulent eddies is known to be better represented by Levy-type distributions 
than by Gaussian/Maxwellian distributions. The collisional relative velocities p.d.f.s could therefore 



differ significantly from the one given in equation ( 26 ) . Levy- type distributions are defined by much 



more slowly decaying tails; this is likely to affect the particle size distribution evolution significantly 
(although quantitatively rather than qualitatively). 

Finally, of course, we must remember that the models studied in Section [5] are local models, 
in the sense that they neglect the radial flux of particles in and out of the domain considered. This 
approximation is adequate a long as the mass flux in and out of a mass bin caused by radial drift is 
small compared with that caused by coagulation and/or fragmentation. As we saw, including mass 
transfer and velocity p.d.fs is always sufficient, in a local model, to trigger the growth of a large- 
particle population. In a global disk model, however, this effect will have to be sufficiently strong 
to overcome the effect of radial drift on the large particle growth. Interestingly, however, we now 
understand what controls the growth rate beyond the fragmentation barrier, so the observations 
of the ubiquity of exoplanets, as well as of the presence of a cm-size particle population at large 
radii, are still reconcilable with our type of model, and will help place constraints on assumed 
sticking and fragmentation parametrization and on the particle velocity p.d.f. selected. Finally, 
note that taking into account the effect of radial drift in our velocity p.d.f.s enhanced particle growth 
very significantly compared with the similarly local Maxwellian-only type of model proposed by 



Windmark et al. (2012b). As such, it is much more likely to yield growth to large sizes despite drift 



than theirs, and should therefore be used preferentially. 



6.2. Observational perspective 

Various observational studies of protoplanetary disks have inferred the presence of a wide 



distribution of particle sizes from submicron-sized dust grains to cm-sized pebbles (Weintraub 



eTaLl[T989l |Dutrey et aE1[T996l |Holland et aTp998| |Beckwith|[l999l |Lada et al.|[2000| |Wilner et al 



2005 Sicilia-Aguilar et al. 2006 Andrews & Williams 2008; Lommen et al. 2009). Intriguingly, 



no conclusive evidence has been found for a correlation between the age of primordial disks and 
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the properties of their population of small dust grains (as revealed by silicate features: Kessler- 
Silacci et aL|2006[ |Furlan et aLl|2009[ |01iveira et "al]|2010[ ). |Williams & Cieza| ( |2011[ ) attribute this 
to a balance between grain growth and destruction on the one hand, and between crystallization 
and amorphization on the other, concluding that this balance seems to persist throughout the 
duration of the primordial disk stage (at least in the surface layers of the disk). The consequence 
for the scenario in which planet formation proceeds via dust growth is the coexistence of particles 
growing to planetary sizes with a population of small grains that are resupplied by continuous 
fragmentation. Future observations of planets embedded in a young protoplanetary disk or the 



confirmation of candidate planets in transition disks, such as T Cha (Huelamo et al. 2011) or 



LkCal5 (Kraus & Ireland 2012), would support this picture. 



In the meantime, other proxies have been sought for evidence of ongoing planetary formation. 
Unfortunately, the failure of theoretical collisional models to produce results in which large and 
small particles co-exist have, until now, hindered these efforts. The new model presented in this 
paper paves the way to resolving this problem. Our physically motivated coagulation and fragmen- 
tation kernels are able to simultaneously establish a sustained growth of particles and to preserve a 
micron-sized dust population in a protoplanetary disk. The particle size distribution that naturally 
emerges is one in which the small and large particle populations have notably different power-laws: 
one that is dominated by fragmentation, and one that is dominated by coagulation via sweeping. 
Since the upper size-cutoffs and total mass in each population are strongly dependent on the model 
considered, observations may help to rule out certain aspects of the parameter space that affect the 
theoretical results, enabling the model to be constrained. 

From a qualitative point of view, the existence of two particle populations far from the central 



star agrees well with detailed modeling and cm-observations of the TW Hydrae disk by |Calvet et al. 



(2002) and Wilner et al. (2005). Furthermore, we expect large and small particles to coexist at 



different radial positions in the disk, albeit with a strong variation in the maximum size achievable 
by the larger particles - in other words, one may anticipate strong radial gradients in inferred 



particle growth beyond mm-size. This result is interesting in the light of the work of Guilloteau 



et al. (2011), who reported observational evidence for a radial dependence of the grain size in 



protoplanetary disks. We expect that more detailed future observations probing different parts of a 
disk will be able to determine whether or not two particle populations do indeed coexist at multiple 
radii in a disk, as suggested by our model. In any case, in light of our results, observers should 
always consider the possibility that two populations of particles exist when fitting spectral energy 
distributions of disks, rather than a single continuous population from small to large sizes. 
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A. Appendix A: Derivation of the relative velocity distribution function 



To derive equation (26) from (25), we expand A^ in spherical coordinates (to be specified, see 
below) as 



A 



raj 



A 



Aij sin 9 cos < 
Aij sin 9 sin <; 
cos 9 , 



(Al) 



and integrate equation (25) over a spherical surface of radius A^-: 

"7T r-lix 



/*7T pZTT 

p(Aij)dAij = / / p 3£) (Ay)A?sin 
io Jo 



(A2) 



Assuming isotropy in the particle dispersion (see Section 3.2), one can re-write p^£>{Aij) as 



p 3D (Aij) 



(2^)3/2^3 



cxp 



A?, + (A 



D\2 
ij) 



2 4 



oh 



(A3) 



where Af? 



(It,; - Uj,Vi - Vj,Wi - W 



•j) and Af- = |A^| (see equation (18)). The integral over the 
spherical surface of p 3 £>{Aij), written in this form, is much simpler if the expansion in spherical 



This preprint was prepared with the AAS IATj^X macros v5.2. 
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coordinates ( Al ) is cleverly chosen with Af- defining the polar axis, so that Ay • A® = Ay A^ cos 9 
We then have 

p(A ij )dA ij = — NQ/0 " q exp j ^ ,, - - | / cxp 



(2vr) 3 / 2 4 



A^ + (Ag) 2 ' 



'A 4j Agcos(9 N 



sin 8 d0 , 



(A4) 



which can be integrated to recover (26). 



